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We study the Lyapunov instability of a two-dimensional fluid composed 
of rigid diatomic molecules, with two interaction sites each, and interacting 
with a WCA site-site potential. We compute full spectra of Lyapunov ex- 
ponents for such a molecular system. These exponents characterize the rate 
at which neighboring trajectories diverge or converge exponentially in phase 
space. Qualitative different degrees of freedom - such as rotation and trans- 
lation - affect the Lyapunov spectrum differently. We study this phenomenon 
by systematically varying the molecular shape and the density. We define and 
evaluate "rotation numbers" measuring the time averaged modulus of the an- 
gular velocities for vectors connecting perturbed satellite trajectories with an 
unperturbed reference trajectory in phase space. For reasons of comparison, 
various time correlation functions for translation and rotation are computed. 
The relative dynamics of perturbed trajectories is also studied in certain sub- 
spaces of the phase space associated with center-of-mass and orientational 
molecular motion. 

PACS numbers: 05.45. -hb, 02.70.Ns, 05.20.-y, 66.90.-hr 



I. INTRODUCTION 

The chaotic molecular motion in fluids and (nonlinear) solids has been mainly studied 
in the past in terms of correlation functions or related power spectra of assorted dynamical 
variables. There is, however, a more fundamental point of view: The basic underlying 
dynamical processes are collisions (interactions) between particles with convex potential 
surfaces. As a consequence, the phase space trajectory is highly unstable which is reflected 
in a very sensitive dependence on initial conditions. This phenomenon is characterized in 
terms of the set of Lyapunov exponents {A/}, / = 1, ■ ■ ■ , L, usually ordered from the largest 
to the smallest. The largest exponent Ai describes the time-averaged logarithmic rate at 
which nearby phase-space trajectories separate. The sums of exponents Z]i=i Aj describe the 
expansion or contraction rates of /—dimensional phase-space objects. Thus, the Lyapunov 
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exponents represent the time-averaged local deformation rates in the neighborhood of a 
phase-space trajectory specified by the time evolution of specially-selected perturbation 
vectors 6i{t). From a practical point of view the precise orientation of the set of initial vectors 
{5;(0)} is not known, nor is it needed for the determination of the A;. This is discussed in 
more detail in Section II. The total number L of exponents is equal to the dimensionality of 
the phase space, and the whole set of exponents is referred to as the Lyapunov spectrum. 
For the evaluation of all L exponents the simultaneous integration of L x (L + 1) first-order 
differential equations is required. The method has therefore been restricted to rather low- 
dimensional dynamical systems in the past. The feasibility of such studies for many-body 
systems has been demonstrated by Hoover and Posch who investigated fluid systems 
with up to 100 atoms in two dimensions. 

Lyapunov spectra of Hamiltonian systems in thermal equilibrium exhibit a pronounced 
symmetry which helps to reduce the number of differential equations to L(L + 2)/2: for each 
positive exponent there exists another negative exponent with equal absolute magnitude. 
This is referred to as Smale pairing |l| or conjugate pairing 0,^. It is a consequence of the 
symplectic nature of the equations of motion |^ , which means that the phase flow - viewed 
as a canonical transformation of the phase space onto itself - leaves the differential two-form 
YliiLx dpi A dqi invariant. Here the sum is over all degrees of freedom, and Pi and denote all 
components of particle momenta and positions, respectively. Due to this pairing symmetry 
the calculation may be restricted to the positive exponents, thus permitting the simulation of 
more complex and larger systems. With the available computer hardware systems with up to 
400 degrees of freedom may be simulated at present by far exceeding the complexity one 
usually encounters when studying dynamical systems with a low-dimensional phase space 

i- 

For atomic fluids one flnds that the shape of the Lyapunov spectrum changes qualita- 
tively if the density is isothermally increased from that of a dense gas to solid densities [§ , 
and that the largest Lyapunov exponent Ai exhibits a maximum at the solid-fluid phase 
transition density 0. From a simulation of thermostatted butane molecules one may 
further deduce that the largest contribution to Ai - and therefore the largest source for chaos 
- is due to the torsional motion around the central CC-bond. The other degrees of freedom 
contribute mainly to the smaller exponents |jlO[|. This conclusion has been reached by ob- 



serving the variations in the Lyapunov spectrum of a Nose-Hoover-thermostatted molecule, 
if various contributions to the molecular Hamiltonian specifying different degrees of freedom 
are selectively switched off. Considering these examples we may expect that also for molecu- 
lar fluids the stretching and contraction properties are very differently affected by translation 
and rotation, and that - as a consequence - the shape of the Lyapunov spectra will strongly 
depend on the state of the system. To clarify this point we report in this paper flrst results 
of a molecular dynamics simulation study of a simple planar molecular fluid in equilibrium 
consisting of N rigid homonuclear diatomic molecules with anisotropy d/a. Here, d is the 
flxed bond length separating the two interaction sites of a molecule, and a is an atomic-size 
parameter d, or the related moment of inertia, and the molecular number density n* are 
two independent parameters with which the emphasis may be conveniently shifted between 
translational and rotational dynamics, at the same time varying the respective relaxation 
times. We measure the associated changes in the (in)stability properties of the phase space 
trajectories. 
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In Section II we define our model fluid and indicate our method of simulating the reference 
trajectory. In Section III we shortly outline the methods for the evaluation of Lyapunov 
spectra. We deflne so-called "rotation spectra" for the vectors 6i connecting neighboring 
trajectories in phase space It turns out that it is useful to study also the projected motion 
of these vectors in subspaces of the phase space such as the respective configurational and 
momentum spaces for translational and rotational motion. To make contact with the more 
traditional analysis of molecular dynamics in fluids, we also evaluate correlation functions 
pertinent for these degrees of freedom. In Section IV we summarize our results. They 
are discussed further in the concluding Section V, where we speculate also about possible 
interpretations of the Lyapunov spectra in terms of suitably deflned collective modes in 
dense fluids. 



II. THE MODEL FLUID 

Simulations were performed for a purely classical system consisting of = 18 rigid 
homonuclear diatomic molecules in a two-dimensional square box of area (volume) V and 
with periodic boundary conditions. The two interaction sites of each molecule are separated 
by a rigid distance d along the molecular axis, and sites on different molecules separated by 
r interact with a purely repulsive Weeks-Chandler-Anderson potential. 
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, r > 2i/V • ^^"^ 

The total intermolecular potential energy $ is taken to be pairwise additive. Reduced units 
are used throughout, for which e, a and the atomic mass m are unity. The equations of 
motion for the 2N interaction sites were augmented with constraint forces keeping the bond 
length d for each molecule flxed |Tl| ]: 

Pi = fi + /i(q2-qi) (2) 

P2 = f2 - /u(q2 - qi). 

Here, q^, are the respective position and momentum vectors of the two interaction sites 
of a molecule, a = 1, 2, and = — Va$ is the force on site a. Furthermore, 

_ (P2 - PiY/m + (q2 - qi) ■ (f2 - fi) , . 

2(q2 - qO^ 

is the Lagrange multiplier for the holonomic constraint (q2 — qi)^ — (i^ = for this molecule. 
The initial conditions were chosen such that the total center of mass velocity was zero. 

For the evaluation of the Lyapunov exponents and the ensuing discussion it is use- 
ful to represent the state of the system by the 6A^-dimensional state vector T{t) = 
{xi,yi,Xi,yi,rii,fii}, i = l,N, where Xi,yi and Xi,yi are the center of mass coordinates 
and momenta, respectively, of molecule i, and rji is the angle the molecular axis of i makes 
with some arbitrary direction. i)i is the corresponding angular velocity. For each molecule i 
these quantities may be easily computed from the instantaneous positions and momenta of 
the interaction sites. 
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III. LYAPUNOV EXPONENTS AND ROTATION NUMBERS 



The equation of motion for the state vector r{t) is conveniently written as an autonomous 
system of first-order differential equations: 

r(t) = G(r(t)). (4) 

Its solution defines a flow r{t) = $j(r(0)) in phase space. Let r(0) be the initial condition 
of a reference trajectory, Ts{0) the initial point of a neighboring perturbed trajectory, and let 
these two points be connected by a parametrized path Co(s) with a perturbation parameter 
s such that lim5_»o ^s{0) = r(0). At time t these points will be mapped by the flow into the 
points T{t) = $t(r(0)) and r,(t) = <l>j(r,(0)), and the path Co(s) into Cj(s). Now we can 
define a finite-length tangent vector at t = 0, 

.(0) = ii„r.(0)-r(0) 

associated with an initial perturbation Ts{0) — r(0) of the reference trajectory in phase 
space. As time goes on, this perturbation develops into Ts{t) — T{t), and the associated 
tangent vector becomes 

.W^UojMl^. (6) 

s~*0 s 

The change of the length of this vector during the time interval t determines the stability of 
the reference trajectory due to the initial infinitesimal perturbation. d{t) may be viewed as 
a vector comoving and corotating with the phase flow in the immediate neighborhood of the 
phase point. It specifies a direction in phase space which varies with time. The equations 
of motion for d{t) are obtained by linearizing the original motion equations (^), 

(5(t) = D(r(t))-5(t), (7) 

where D(r) = [{d/dr)G{T)] is an L x L matrix and is referred to as the stability matrix. 
The formal solution of may be written as 

S{t) = M{t;5{0))-5{0), (8) 

where the operator M is a time ordered exponential of J T){{r{t'))dt' . According to the 
multiplicative ergodic theorem by Oseledec 0,0] the Lyapunov exponents 

Xi = lim -ln|M(t;5(0))uz| 

= Inn iln|ul-H(t;5(0))-u;| (9) 

exist for mixing systems and are independent of the initial conditions. In this equation 
u; denote the L orthonormal eigenvectors of the real and symmetric matrix H(t; 6(0)) = 
M(t; (5(0))^M(t; (5(0)), where f means transpose. They may be taken as the basis for an 
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arbitrary inital vector 6(0), for which the long-time behavior is ultimately determined by 
that vector component d{0) ■ u; with the largest associated Lyapunov exponent. 

The classical algorithm for the calculation of the complete spectrum of Lyapunov ex- 
ponents due to Benettin et. al. [|l5|,|l6l and others requires the simultaneous 



integration of the original reference system and of L sets of the linearized equations (|^ 
for L initial tangent vectors Si{0) taken to be orthonormal. However, due to the stretch- 
ing and folding operations of the phase flow, these vectors will not stay orthonormal for 
t > but will be stretched and rotatad into the direction of the largest phase-space ex- 
pansion corresponding to the largest exponent, and eventually diverge. This is prevented 
by reorthonormalizing the vectors after every few time steps. The Lyapunov exponents are 
determined from the time-averaged contraction/expansion factors for the vector norms. 
A conceptual refinement of this algorithm has been first proposed by Hoover et. al. [p!2|,p|. 



and independently by Goldhirsch et. al. [jT9l. It has again been reinvented since then pO[ . 



In this method the vectors 6i are constrained to remain orthonormal for all times t > by 
the introduction of constraining forces added to the right-hand side of the linearized motion 
equations: 

di{t) = B{T{t))di-j2\viSu (10) 

i'=i 

The Xi'^i are time dependent Lagrange multipliers and are determined from the orthonor- 
mality conditions (<5; -61') = dw: 

The Lyapunov exponents are given by the time-averaged diagonal multipliers 

Xi = {\n{t)). (12) 

This algorithm is equivalent to the classical algorithm with continuous reorthonormalization 
and yields solutions for the tangent vectors 5i{t) identical to the classical ones, if in the latter 
case the vectors are reorthonormalized after every time step. However, the equations ( pAj ) 
together with (arbitrary) initial conditions may be considered as the defining equations for 
the orthonormal vectors 6i{t). They emphasize another property of tangent-space dynamics 
which has been virtually ignored up to now |T9|J^, namely the rotation of the orthonormal 



tangent vectors. 

We have seen that the solution of Equ. (0) constitutes an orthonormal set of vectors 
which continuously change their orientation in tangent space with instantaneous angular 
velocities The vectors are forced to remain orthogonal to each other through the 

force terms proportional to the off-diagonal Lagrange multipliers in equation (P). If viewed 
in phase space, these orthonormal vectors move with the state point along the reference 
trajectory and simultaneously reorient such that 61 always turns toward the direction of 
fastest phase-space expansion, ^2 into a perpendicular direction with second largest growth, 
and so forth. As a measure for this unitary rotation we have computed an averaged angular 
velocity for each vector 5i defined by 
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1 ^ cos-\di{Q-di{t^ + At)) 



_ 1 ^ |A9Ktn)| 

where AG/(t„) is the angle by which the unit vector 5i reorients during a time step At at 
time t„ = riAt, and Nt is the number of time steps of the simulation. We refer to these 
numbers as "rotation numbers", and to their whole set as the "rotation spectrum" [P].Thus, 
uji is the time-averaged modulus of the angular velocity for the reorienting vector 6i{t) in 
phase space. 

We stress that - unlike the Lyapunov exponents - the rotation numbers defined in (|T3|) 
depend on the metric of the phase space and of the coordinate system used. There is no 
multiplicative ergodic theorem for these quantities, although they are independent of the 
choice of the initial conditions. In this respect they are on the same level of theoretical 
significance as the instantaneous finite-time Lyapunov exponents, which also depend on the 
choice of the coordinate system |0. This issue clearly needs further investigation. In spite of 
these theoretical restrictions, the rotation numbers still convey important information about 
the phase space dynamics. For example, for isothermal scans through order-disorder phase 
transitions the rotation numbers increase monotonously with density, whereas the maximum 
Lyapunov exponent exhibits a pronounced maximum at the transition density [^0. The 
rotation numbers are also largest for indices belonging to the smallest exponents believed to 
be associated with the stability of collective modes in a fluid. We therefore speculate that 
the rotation spectra may prove useful for establishing a link between the linearized dynamics 
in tangent space and more traditional descriptions of systems in terms of (collective) modes. 
Examples for rotation spectra are given in Section IV. Related numbers have been defined 
by Ruelle |2l[] and applied to two-dimensional maps by Lambert et. al. |p2 |. 



It is instructive to view the dynamics of the (5-vectors not only in the whole phase 
space but also in certain qualitatively different subspaces associated with special degrees of 
freedom. Since the phase space is a product of the center-of-mass configuration space Q, of 
the respective momentum space P, of the angular orientation space VL for the molecular axes, 
and of the associated angular momentum space Po, also the tangent space is decomposed 
into respective subspaces TQ, TP, Tf2, and TPq. We consider projections of the 5i vectors 
onto TX e {TQ, TP, TO, TPn}: 

5x,i = V{TX)6i. (14) 

The projection operator V(TX) may be represented as a diagonal matrix with elements 
VaaiTX) equal to unity, if the a-axis of 8i belongs to TX, and zero otherwise. We compute 
the time-averaged squared lengths of these projected vectors: 

4,^ = {Sx,i ■ Sx,i) (15) 

and refer to them as "mean-squared X-components" of Si. Of course, for each / they add up 
to unity if summed over TX. They are a measure for the probability of a vector di of pointing 
into a direction of tangent space belonging to the subspace TX. They turn out to be very 
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helpful for the interpretation of our results. A related quantity, cos^ a*^"^) = been 
discussed recently by D'Alessendro and Tenenbaum who refer to a^^^ as a coherence 
angle. It represents an effective angle between the subspace TX and the maximum-expansion 
subspace. 

We have mentioned already that the constrained orthonormal- vector method of equations 
( p!0| - |T^) is conceptually very useful. However, any numerical solution of these equations 
involves quite extensive vector- matrix operations. Since the tangent vectors Si{t) obtained 
in this way are identical to the vectors generated by the classical method with continueous 
reorthonormalization, equations (|10|-|12]) do not offer an improvement from a numerical point 
of view. For our numerical work we have therefore used the classical method of Benettin 
in combination with Gram-Schmidt reorthonormalization after every time step. To avoid 
the algebraic complexity we employed a finite-difference variant of the classical algorithm: 
The tangent vectors in and (P) were approximated by finite offset vectors between two 
phase-space trajectories corresponding to a finite s = 0.0001. As many sets of the original 
motion equations as the required number of Lyapunov exponents were integrated for 
the determination of the satellite trajectories, and Gram-Schmidt reorthonormalization was 
carried out after every time step. A fourth order Runge-Kutta algorithm with a reduced 
time step At = 0.001 was used for the integration. Every few time steps the molecular 
bonds were rescaled to their precise length to compensate for numerical inaccuracies. In all 
runs the trajectories were followed for at least 400 time units. 



IV. RESULTS 

The reduced molecular number density n* = Na'^/V was chosen large enough to ensure 
sufficient coupling between translational and rotational degrees of freedom. Two series of 
simulation runs with n* equal to 0.4 and 0.5 were performed. For a given n* the molecu- 
lar shape was varied by considering systems with different molecular bond length d. The 
anisotropy parameter d/a was varied between 0.2 and 1.0, and the corresponding variation 
of the Lyapunov spectrum was determined. Since n* does not account for the molecular 
anisotropy, we characterize our model systems also by an anisotropy-dependent density pa- 
rameter n'^ = Na{a + d)/V which is, roughly speaking, the ratio of the occupied volume to 
the total, n'^ becomes equal to n* for isotropic particles. Some of our results for the maxi- 
mum Lyapunov exponent Ai and for the Kolmogorov entropy are summarized in Table 

L /2 

1 for n* = 0.4, and in Table 2 for n* = 0.5. Here, the Kolmogorov entropy = J2i=i is 
given by the sum of the non-negative exponents. For all these simulations the kinetic energy 
per molecule was equal to 2, 2/3 for each translational and rotational degree of freedom. 
We monitored this partition of the kinetic energy among the translational and rotational 
degrees of freedom throughout the simulation. We found that equipartition among these 
degrees of freedom holds to within 3% for the whole range of densities, bond lengths and 
temperatures. 

The very different n'^-dependence of Ai from that of hk already suggests that the Lya- 
punov spectrum is considerably affected by the bond length. This may be verified from 
Figure 1 in which all spectra with identical n* = 0.4 but different d are displayed, and from 
the analogous Figure 2 for n* = 0.5. Of course, each spectrum consists of discrete points 
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only, which arc located at the nodes of the connecting lines. 

The index / along the abscissa merely numbers the exponents (or degrees of freedom), 
with / = 1 for the maximum exponent, 51 for the smallest positive , and 108 for the most 
negative exponent. For the construction of Figs. 1 and 2 only the positive branches (54 
exponents) of the Lyapunov spectra were calculated and displayed. Due to the Smale-pairing 
symmetry for symplectic systems mentioned in the Introduction, the negative branch is 
obtained by reversing the sign of the positive branch. 6 of the exponents must vanish due to 
the 5 constants of the motion - energy, center of mass, and linear momentum - and the fact 
that any perturbation in the direction of the phase flow adds another vanishing exponent. 
Vanishing exponents and their associated phase space directions have indices 52 < Z < 57. 
However, for short bond lengths the exponent A52 in Figs. 1 and 2 does not vanish exactly. 
This is an undesirable consequence of the periodic rescaling of the molecular bonds, but 
does not affect the shape of the spectra. 

For a few selected systems containing molecules with the longest {d = a) and the shortest 
bonds {d — 0.2a) we computed also the full spectrum of 108 exponents together with the 
respective rotation numbers and mean-squared X-components d\ ^. For the longest molecule, 
d = a, full rotation spectra are depicted in Figure 3 for the densities n* = 0.4 and 0.5. Mean- 
squared X-components are displayed in Fig. 4 for n* = 0.5, d = a, and in Fig. 5 for n* = 0.5, 
d — 0.2(7. These results will be discussed further in the following Section. 

In order to relate these findings to more traditional descriptions of the dynamics, we 
calculated also various time correlation functions referring to translational and rotational 
motion of our model system. The normalized velocity autocorrelation functions for molecular 
center-of-mass motion, Cvv{t) = (v(t) • v(0)) / (v(0) ■ v(0)), for the systems with the larger 
density n* = 0.5 are shown in Fig. 6. Their time integral is related to the translational 
diffusion coefficient. Furthermore, the reorientational dynamics of linear molecules is usually 
discussed in terms of the correlation functions for the spherical harmonics of rank / of the 
angles specifying the orientation of the molecular axis. In three dimensions some of these 
correlation functions may be accessible by experiment, such as infrared absorption (/ = 1), 
Raman scattering (/ = 2), and neutron scattering (a weighted sum of / > terms). For our 
planar model these rotational correlation functions are written as 

Ci{t) = {Pi{cos{Qm, (16) 

where 0{t) — r](t) — r){0) is the angle the axis of a molecule reorients during the time t, and 
Pi denotes a Legendre polynomial of rank /. In Figure 7 the functions Ci{t) = (cos(0(t))) 
are depicted for the case n* = 0.5 and the whole range of bond lengths. Similar results have 
been obtained also for C2 but are not shown here. 

V. DISCUSSION 

The computation of L Lyapunov exponents requires the simultaneous integration of 
L(L + 1) first-order differential equations. Thus, the number of molecules of a system 
accessible to present-day computation is much smaller than one is used to for traditional 
simulations of dynamical properties such as correlation functions or transport coefficients. 
For that reason our system consists of only N — IS two-center molecules in a plane, and the 
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simulation box has a typical side length of Qa (for n* = 0.5). For such small systems it is 
difficult to distinguish between fluid and solid phases, although already systems of only two 
disks exhibit a first-order phase transition Nevertheless, we shall use these expressions 



to convey the general idea. 

Before entering the discussion we stress again that 6 of the Lyapunov exponents vanish 
for reasons given above. They are located in the middle of the spectrum, 52 < / < 57, 
and the associated phase-space directions belong to the central manifold. Since Gram- 
Schmidt reorthonormalization has no ordering effect on the directions of their <5;-vectors the 
respective rotation numbers in Fig. 3 and the squared X-components in Figs. 4 and 5 have 
no meaning for 52 < / < 57. We have nevertheless included these points in the Figures to 
draw attention to this peculiarity. 

An inspection of the center-of-mass velocity autocorrelation functions in Figure 6 reveals 
that backscattering - typical for a dense fiuid or solid - occurs only for the largest bond 
length. For this system Ci in Fig. 7 does not decay, and the orientation of the molecules 
persists for a long time. We conclude that this state is a solid, whereas all other systems are 
fluids. This is also reflected in the shape of the Lyapunov spectra, which for d = 1.0 and 
n* = 0.5 in Fig. 2 is very similar to that of a two-dimensional solid formed with isotropic 
particles B. Also the empirical power law 



A,/Ai 



(|- 2 -0/(^-3) 



n/3 



(17) 



for l<Z<L/2 — 2 works reasonably well with an exponent /5 = 3/2 0. Ai is the maximum 
exponent. As expected, the shape of all other spectra conforms closely to that of atomic 
fluids [0], but a representation in terms of such a power law is not appropriate. 

The maximum Lyapunov exponent Ai describes the time evolution of the small pertur- 
bation which grows fastest in phase space. It is a local quantity in the sense that it depends 
basically on the fastest dynamical events taking place in the system, namely collisions, for 
which the velocities and angular velocities change sign. It is a reasonable assumption that the 
fastest phase-space growth takes place in the linear and angular momentum subspaces. The 
mean-squared X-components introduced in Section 3 support this view. Simulations 
for atomic systems have revealed that <5i is located almost fully in the momentum related 
subspace TP of tangent space. For the linear molecules studied here the angular-momentum 
related subspace TPq of tangent space turns out to be the most important. Fig. 4 shows 
that on the average 69% of the squared length of <5i for the elongated molecular case {d = 1) 
is contributed by TPq. This number even rises to 96% for the more freely rotating case in 
Fig. 5 {d = 0.2). We conclude that the main reason for the instability of the phase space 
trajectory is due to the anisotropy of the pair potential, and it is mostly accumulated in the 
angular momentum subspace. 

The situation changes completely when we consider ^(-vectors for / > 1 pointing into 
less-violently expanding or even compressing phase-space directions. Figs. 4 and 5 reveal 
that the linear-momentum subspace becomes dominant for, say, 20 < / < 51, still associated 
with positive exponents. In this range the prominent contributions to the Lyapunov spectra 
comes from translational modes which one is tempted to associate with generalized "hydro- 
dynamical" modes with small but flnite wave vectors. An approximate representation of 
these translation-dominated exponents in terms of the power law (|T^ leads to exponents 
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(3 < 1 similar to that of atomic fluids 0. It follows that the positive curvature of the Lya- 
punov spectra in Figs. 1 and 2 for the most positive exponents (small /) is essentially due to 
contributions from rotational degrees of freedom which make themselves felt more distinctly 
for small d associated also with small moments of inertia. 

It is interesting to note that the stable phase-space directions corresponding to the neg- 
ative Lyapunov exponents are dominated by the configurational subspace Q and, to a lesser 
extent, by the orientation-angle subspace Q (Fig. 4). The significance of this is not clear to 
us. 

That the maximum exponents of Figs. 1 and 2 increase with density n* for fixed d and 
constant temperature is easily explained by the increase of the collision frequency. The 
relative maximum exhibited by the positive exponents of Fig. 2 as a function of the bond 
length is less obvious, since an increase of d for constant n* also increases the collision rate. 
This maximum is a consequence of the phase transition eventually leading to a solid for 
d = 1. In the less-dense case of Fig. 1 this maximum is expected to occur for bond lengths 
slightly larger than a. Similar maxima for Ai as a function of density have been observed 
previously for fluid-to-solid phase transitions in atomic-particle systems p]. 

The rotation numbers measure the average speed of rotation of the ^(-vectors in phase 
space. From simple arguments involving the time-reversal invariance of the original motion 
equations (^) and of their linearized version (|^) we expect that the rotation spectra for 
symplectic systems are symmetric such that uJi = uj^^i^i. This symmetry was observed for 
atomic systems §], and is also apparent in Fig. 3 for the linear-molecule case with d = 1. 
The theoretical significance of these numbers is still controversial. Also the A^-dependence 
of these numbers needs to be investigated. 

The Kolmogorov entropy hx is a global measure of the rate with which information is 
generated by the dynamics and, hence, of the disorder in such an equilibrium system. Our 
numerical results tabulated in Tables 1 and 2 confirm the general picture outlined above. 
For n* = 0.5 this parameter varies only very little with the molecular anisotropy as long 
as the system remains fluid (Table 2). It starts to decrease, when the phase transition is 
approached, and becomes significantly smaller for the solid. This transition takes place near 
an anisotropy-dependent density n'^ = 0.8. For the lower density n* = 0.4 a similar transition 
occurs near n"^ = 0.8, although d must be increased beyond 1 for a solid to be observed. 
This "transition density" agrees with the phase-transition density of an isotropic-particle 
system, at which the maximum Lyapunov exponent reaches a maximum 

VI. CONCLUSIONS 

This work is a first attempt to survey the (in)stability properties of phase-space trajec- 
tories for systems of anisotropic molecules. A detailed analysis of the Lyapunov spectra as 
a function of density and molecular anisotropy makes it possible to distinguish - at least 
qualitatively - the contributions from the center-of-mass motion and of molecular reorienta- 
tion. We find that the major contributions to the instability of the phase space trajectory 
come from the rotational degrees of freedom and, in particular, from the angular momen- 
tum variables. Translational center-of-mass motion is much less destabilizing. Although, 
for practical reasons, the systems contain only a few molecules, the influence of the fluid- 
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solid phase transition on the Lyapunov instabihty and the Kolmogorov entropy is clearly 
seen. One might speculate that an even more detailed analysis of the dynamics of the 
tangent vectors in the respective subspaces spanned by the center-of-mass coordinates, the 
corresponding momenta, the orientation angles, and the angular velocities may lead to an 
interpretation in terms of collective modes in many body systems. We are still far away from 
this goal, but it is hoped that recent progress in the methodology of computing Lyapunov 
spectra for systems of hard core particles P5|J2B| J5[] will also be useful for the understanding 
of anisotropic molecular systems. 



VII. ACKNOWLEDGEMENTS 

We thank Ch. Dellago and W.G. Hoover for many stimulating and helpful discussions 
on this and related subjects. This work was supported, in part, by a grant to I. Borzsak 
by the Soros Foundation (S-2374/93), in part, by the Austrian Fonds zur Forderung der 
wissenschaftlichen Forschung, Grant P9677, and, in part, by the Hungarian OTKA, Grant 
F7218. We are also grateful to the staff of the Computer Center of the University of Vienna 
for the generous allocation of computer resources. 



11 



H. A. Posch and W. G. Hoover, Phys. Rev. A 38,473 (1988). 

H. A. Posch and W. G. Hoover, Phys. Rev. A 39, 2175 (1989). 

D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Phys. Rev. A 42,5990 (1990 ). 

S. Sarman, D. J. Evans, and G. P. Morriss, Phys. Rev A 45, 2233 (1992). 

V. I. Arnold, Mathematical Methods of Classical Mechanics, Springer, New York (1983). 

H. A. Posch and W. G. Hoover, Phys. Rev. E 49, 1913 (1994). 

W. G. Hoover, C. G. Hoover, and H. A. Posch, Phys. Rev. A 41, 2999 (1990). 

H. A. Posch, W. G. Hoover, and B. L. Holian, Ber. Bunscngcs. Phys. Chcm. 94, 250 (1990). 

Ch. Dellago, H. A. Posch, and W. G. Hoover, Phys. Rev. E , in print (manuscript EG5641). 

H. A. Posch and S. Toxvaerd, unpubUshed. 

D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids. Academic 
Press, 1990. 

W. G. Hoover and H. A. Posch, Phys. Lett. A, 113, 82, (1985). 

V.I. Oseledec, Moscow Math. Soc, 19, 197 (1968). 

J. P. Eckmann and D. Ruelle, Rev. Mod. Phys. 57, 617 (1985). 

G. Benettin, L. Galgani, A. Giorgilh, and J. M. Strelcyn, C. R. Acad. Sci. (Paris) 286, 431, 
(1978). 

G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, Meccanica 15, 9 (1980). 

I. Shimada and T. Nagashima, Progr. Theor. Phys. (Japan) 61, 1605 (1979). 

A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, Physica D 16, 285 (1984). 
I. Goldhirsch, P. L. Sulem, and S. A. Orszag, Physica D 27, 311 (1987). 
W. E. Wicsel, Phys. Rev. E 47, 3686 (1993) 
D. Ruelle, Ann. Inst. Poincare 42, 109, (1985). 

A. Lambert, R. Lima, and R. Vilela Mendes, Physica D 34, 366 (1989). 
M. D'Alessandro, and A. Tenenbaum, Phys. Rev. E 52, R2141 (1995). 

B. J.Alder, W.G. Hoover, and T.E.Wainwright, Phys. Rev. Lett. 11, 241 (1963). 
Ch. Dellago, and H. A. Posch, Phys. Rev. E 52, 2401 (1995). 

Ch. Dellago, L. Glatz, and H. A.Posch, Phys. Rev. E 52, 4817 (1995). 



12 



Tables: 



d 


0.2 


0.33 


0.5 


0.66 


1.0 




0.48 


0.532 


0.6 


0.664 


0.8 


Ai 


4.11 


3.74 


3.57 


3.59 


3.47 


hx 


111.0 


114.6 


117.2 


120.0 


111.7 



Table 1: Simulation results for a density n* = 0.4. d is the molecular bond length and is 
given in units of a. n'^, the anisotropy-dependent density, is defined in the text. The 
maximum Lyapunov exponent Ai and the Kolmogorov entropy hx are given in units 
of (e/m(72)V2. 



d 


0.2 


0.33 


0.5 


0.66 


1.0 




0.6 


0.665 


0.75 


0.83 


1.0 


Ai 


5.08 


4.56 


4.36 


4.02 


3.16 


hx 


134.1 


135.9 


133.3 


115.0 


65.2 



Table 2: Simulation results for a density n* = 0.5. d is the molecular bond length and is 
given in units of a. n'^, the anisotropy-dependent density, is defined in the text. The 
maximum Lyapunov exponent Ai and the Kolmogorov entropy hx are given in units 
of (e/ma2)V2. 
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Figure captions: 



Figure 1: Lyapunov spectra for a planar fluid of 18 linear diatomic molecules at a reduced 
density n* — 0.4, for five different bond lengths d/a equal to 0.2, 0.33, 0.5, 0.66, 
and 1.0. Only the positive branch (54 exponents, of which 3 vanish) is shown. The 
Lyapunov exponents are given in units of (e/mcr^)^/^. As usual, the index I merely 
numbers the exponents. A/, such that I — 1 refers to the maximum exponent. 

Figure 2: Lyapunov spectra for a planar fluid of 18 hnear diatomic molecules at a reduced 
density n* = 0.5, for five different bond lengths d/a equal to 0.2, 0.33, 0.5, 0.66, 
and 1.0. Only the positive branch (54 exponents, of which 3 vanish) is shown. The 
Lyapunov exponents are given in units of (e/mcr^)^/^. As usual, the index / merely 
numbers the exponents. A;, such that I = 1 refers to the maximum exponent. 

Figure 3: Full rotation spectra (108 rotation numbers each) as defined in Section 3, for 
a planar fiuid of 18 linear diatomic molecules with a bond length d/a = 1, for the 
densities n* equal to 0.4 (plus signs), and 0.5 (diamonds). The rotation numbers are 
given in units of (e/mcr^)^/^. The index I merely numbers the spectral points such that 
1 corresponds to the maximum Lyapunov exponent, 108 to the minimum exponent. 

Figure 4: Mean squared X-components 5x i as a function of the Lyapunov index I for a 
planar fluid of 18 linear diatomic molecules with a density n* — 0.4 and bond length 
d/a — 1. 1 = 1 corresponds to the maximum Lyapunov exponent. The subspaces 
X are the center-of-mass configurational subspace Q (diamonds), the center-of-mass 
momentum subspace P (plus signs), the orientation-angle subspace fl (squares), and 
the angular velocity subspace (crosses). 

Figure 5: Mean squared X-components as a function of the Lyapunov index / for a 
planar fiuid of 18 linear diatomic molecules with a density n* = 0.4 and bond length 
d/a = 0.2. I = 1 corresponds to the maximum Lyapunov exponent. The subspaces 
X are the center-of-mass configurational subspace Q (diamonds), the center-of-mass 
momentum subspace P (plus signs), the orientation-angle subspace Q (squares), and 
the angular velocity subspace Pa (crosses). 

Figure 6: Normalized velocity autocorrelation functions for a system of 18 linear diatomic 
molecules with a density n* = 0.5 and different bond lengths d/a varying between 0.2 
and 1 as indicated. The time is given in units of {ma^y^^. 

Figure 7: Orientational correlation functions Ci{t) = (Pi(cos 0(t))) for a system of 18 
linear diatomic molecules with a density n* = 0.5 and different bond lengths d/a 
varying between 0.2 and 1 as indicated. The time is given in units of (mcr^)^/^. 
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